clear all
addpath('../functions')
% before runinng. Run stata program to generate csv. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Main Results
for s = 1:500
    filename = sprintf('../data/bootstrap_samples/csv/UK_Budget_boot_%d.csv', s);
    savename = sprintf('../data/bootstrap_samples/boot/UK_Index_%d', s);
 
    run("../subroutines/MatlabDataConstruction_boot.m")
    [U_vec] = CalMoneyMetric(I_vec, B_vec, pvec,1);
    save(savename,'U_vec','I_vec');
end
%%
load ../data/DataforMatlab/IP.csv
load ../data/DataforMatlab/torn.csv
ns = 500;
    filename = sprintf('../data/bootstrap_samples/boot/UK_Index_%d', 0);
    load(filename)
    N = size(I_vec,2);T=size(I_vec,3);
    run("../subroutines/sub_make_percentile.m")
    OuryaxisUK_base = OuryaxisUK;  
    OurxaxisUK_base = OurxaxisUK;
    EVdivisia_base  = RC; 

OuryaxisUK_vec = zeros(size(OurxaxisUK_base,1),1,ns);

for s = 1:ns
    filename = sprintf('../data/bootstrap_samples/boot/UK_Index_%d', s);
    load(filename)
    F_X = makefunction(log(I_vec(:,:,end)*52),log(U_vec(:,:,end)*52));
    OuryaxisUK_vec(:,:,s) = exp(F_X(log(OurxaxisUK))); % plot at the same points
end

OuryaxisUK_UB = prctile(OuryaxisUK_vec,97.5,3);
OuryaxisUK_LB = prctile(OuryaxisUK_vec,2.5,3);

%%Figure5

axoptions={'scaled ticks = false',...
           'y tick label style={/pgf/number format/.cd, fixed, fixed zerofill,precision=0, set thousands separator={}}',...
           'x tick label style={/pgf/number format/.cd,precision=1, set thousands separator={}}',...
           'legend style={font=\normalsize}'}; %\scriptsize

figure('name','Figure 5','NumberTitle','off')
h = loglog(OurxaxisUK_base ,OuryaxisUK_base,'Color', [0, 0.4470, 0.7410],'LineWidth', 2); hold on

err = [OuryaxisUK_UB';OuryaxisUK_LB'];
hold on
patch([OurxaxisUK_base', fliplr(OurxaxisUK_base')], [err(1,:), fliplr(err(2,:))], 'b', 'FaceAlpha',0.5,'LineWidth', 0.2)   % Shaded Confidence Intervals
loglog(OurxaxisUK_base ,(EVdivisia_base(:,end)),'Color', [0.8500, 0.3250, 0.0980],'LineWidth', 2);
g1 = gca;
xlim([min(OurxaxisUK_base), 100000])
ylim([min(OuryaxisUK_base), max(OuryaxisUK_base)])
xticklabels({'10,000', '100,000'})
yticklabels({'1,000', '10,000'})

ax = ancestor(h, 'axes');
ax.XAxis.Exponent = 0;
ax.XAxis.TickLabelFormat = '%.0f';
xticklabels({'10,000','100,000'})
yticklabels({'1,000','10,000'})
xlabel('Income in 2017')
ylabel('1974 Income that gives the same utility as in 2017')
legend('Money metric','Confidence interval','Real consumption','Real Consumption','Location','northwest')
%matlab2tikz('../fig/Money_boot.tex','extraAxisOptions',axoptions);



% Function to interpolate vectors with NaN values
function F_X = makefunction(X,Y)
    indexToKeep = ~isnan(X) & ~isnan(Y);
    X = X(indexToKeep); Y = Y(indexToKeep);    
    [X, indexUN] = unique(X);
    Y = Y(indexUN);
    F_X = griddedInterpolant(X,Y,'linear','none');
end

